Twin analysis using TwinAnalysis

Ivan Voronin

January 12, 2021

Contents


The R code from this presentation is here: https://raw.githubusercontent.com/IvanVoronin/TwinAnalysis/master/slides/intro/intro.R

Twin method: Fundamentals

ACE

ADE

Univariate ACE

Univariate ACE

Univariate ACE

Univariate ACE

Univariate ACE

TwinAnalysis

Univariate toy data

# Making univariate toy data
library(dplyr,
        quietly = TRUE)

mz_data <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0),
  Sigma = matrix(c(1.0, 0.8,
                   0.8, 1.0),
                 nrow = 2)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'X2'))

dz_data <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0),
  Sigma = matrix(c(1.0, 0.5,
                   0.5, 1.0),
                 nrow = 2)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'X2'))

data <- rbind(data.frame(zyg = 'MZ', mz_data),
              data.frame(zyg = 'DZ', dz_data)) %>%
  mutate(
    zyg = factor(zyg, levels = c('MZ', 'DZ'))
  )

head(data)
##   zyg         X1         X2
## 1  MZ  1.6139277  1.7912405
## 2  MZ -0.7454432 -0.4165514
## 3  MZ -1.2403718 -0.6782260
## 4  MZ -0.9808771 -1.0968942
## 5  MZ -0.3770839 -0.2187076
## 6  MZ -0.6517744 -0.4693972

Descriptives

# Installing the packages
# install.packages('devtools')
# library(devtools)
# install_github('ivanvoronin/mlth.data.frame')
# install_github('ivanvoronin/TwinAnalysis')

# Loading the packages
library(OpenMx)
library(TwinAnalysis)

# Descriptive statistics
get_univ_descriptives(data, zyg = 'zyg', vars = 'X')
## $X
##    Twin1-------------------- Twin2-------------------- Cor--------------------------
##    n   M           SD        n   M            SD       r         lowerCI   upperCI  
## MZ 100 -0.06778361 0.9888608 100 -0.075597225 1.022780 0.7881512 0.7001161 0.8525833
## DZ 100  0.06555748 1.0586309 100 -0.008591523 1.234908 0.5030154 0.3402106 0.6365423

Univariate ACE

# Univariate ACE model
model <- univariate_ace(data, zyg = 'zyg', ph = 'X')

model <- mxRun(model)

summary(model)
## Summary of ACE 
##  
## free parameters:
##            name matrix row col    Estimate  Std.Error A
## 1 ACE.mean[1,1]   mean   1   1 -0.01625766 0.06929624  
## 2    ACE.a[1,1]      a   X   X  0.91070049 0.09748841  
## 3    ACE.c[1,1]      c   X   X  0.37057316 0.22770863  
## 4    ACE.e[1,1]      e   X   X  0.46506707 0.03335832  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:              4                    396
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:               1067.62
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/400
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:        275.620               1075.620
## BIC:      -1030.514               1088.813
##       |  Sample-Size Adjusted
## AIC:                 1075.825
## BIC:                 1076.141
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 04:22:46 
## Wall clock time: 0.07753515 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Output - parameters

# Output tables from the univariate ACE
get_output_tables(model)
## $Variance_components
##        A(%)      C(%)      E(%)
## X 0.7010857 0.1160828 0.1828315
## 
## $Raw_variance
##      Total         A         C         E
## X 1.182987 0.8293754 0.1373245 0.2162874

Output - model fit

# Model fit
model <- ref_models(model)

fit_stats(model)
##   model      base ep minus2LL  df    AIC       BIC
## 1   ACE Saturated  4  1067.62 396 275.62 -1030.514
##         CFI       TLI      RMSEA   diffLL diffdf         p
## 1 0.9929781 0.9976594 0.02696322 6.872418      6 0.3328073

Constrained models

# Constrained models
# This is a list of MxModels
nested_models <- def_nested_models(
  model,
  AE = mxConstraint(c[1, 1] == 0),
  CE = mxConstraint(a[1, 1] == 0),
  E = list(mxConstraint(c[1, 1] == 0),
           mxConstraint(a[1, 1] == 0)),
  run = TRUE
)

fit_stats(model, nested_models = nested_models)
## Warning in computeFitStatistics(likelihood, DoF, chi,
## chiDoF, retval[["numObs"]], : Your model may be mis-
## specified (and fit worse than an independence model),
## or you may be using the wrong independence model, see ?
## mxRefModels
##   model      base ep minus2LL  df      AIC        BIC
## 1   ACE Saturated  4 1067.620 396 275.6200 -1030.5137
## 2    AE       ACE  4 1068.249 397 274.2487 -1035.1833
## 3    CE       ACE  4 1095.945 397 301.9449 -1007.4871
## 4     E       ACE  4 1194.151 398 398.1515  -914.5788
##            CFI       TLI      RMSEA      diffLL diffdf
## 1  0.992978093 0.9976594 0.02696322   6.8724180      6
## 2  0.995966781 0.9988477 0.01891893   0.6286784      1
## 3  0.773045631 0.9351559 0.14191883  28.3249192      1
## 4 -0.009349393 0.7476627 0.27995972 126.5314918      2
##              p
## 1 3.328073e-01
## 2 4.278405e-01
## 3 1.025671e-07
## 4 3.342225e-28

Bivariate toy data

# Making bivariate toy data
mz_data2 <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0, 0, 0),
  Sigma = matrix(
    c(1.0, 0.7, 0.6, 0.5,
      0.7, 1.0, 0.5, 0.8,
      0.6, 0.5, 1.0, 0.7,
      0.5, 0.8, 0.7, 1.0),
    nrow = 4)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'Y1', 'X2', 'Y2'))

dz_data2 <- MASS::mvrnorm(
  n = 100,
  mu = c(0, 0, 0, 0),
  Sigma = matrix(
    c(1.0, 0.7, 0.4, 0.2,
      0.7, 1.0, 0.2, 0.4,
      0.4, 0.2, 1.0, 0.7,
      0.2, 0.4, 0.7, 1.0),
    nrow = 4)) %>%
  as.data.frame() %>%
  setNames(c('X1', 'Y1', 'X2', 'Y2'))

data2 <- rbind(data.frame(zyg = 'MZ', mz_data2),
               data.frame(zyg = 'DZ', dz_data2)) %>%
  mutate(
    zyg = factor(zyg, levels = c('MZ', 'DZ'))
  )

head(data2)
##   zyg         X1         Y1         X2          Y2
## 1  MZ -1.0800889 -1.0275203 -0.1378594 -0.82322342
## 2  MZ  0.8290774  0.9620859  0.9405737  0.88065324
## 3  MZ  0.5680403  0.1000585  1.2833039  0.41037703
## 4  MZ  0.2603521  1.2183755 -0.4669039  0.90408497
## 5  MZ -0.1713043 -0.4440812  0.7967614  0.04784216
## 6  MZ -0.2450042  1.0047931 -0.7089274  0.50509900

Descriptives

# Descriptive statistics
get_univ_descriptives(data2, 
                      zyg = 'zyg', 
                      vars = c('X', 'Y'))
## $X
##    Twin1------------------ Twin2-------------------- Cor--------------------------
##    n   M          SD       n   M           SD        r         lowerCI   upperCI  
## MZ 100 -0.1359043 1.003993 100 -0.06795439 0.9190689 0.5930486 0.4489236 0.7070999
## DZ 100 -0.1636363 1.017908 100 -0.13654309 0.9801215 0.5075549 0.3455896 0.6401541
## 
## $Y
##    Twin1-------------------- Twin2------------------- Cor--------------------------
##    n   M           SD        n   M          SD        r         lowerCI   upperCI  
## MZ 100 -0.02955259 0.9583681 100 -0.1114836 0.9391590 0.7579110 0.6597008 0.8306695
## DZ 100 -0.21072036 1.1037023 100 -0.1693055 0.9790516 0.5489379 0.3951225 0.6728125

Twin correlations

# Cross-twin cross-trait correlations
get_cross_trait_cors(data2,
                     zyg = 'zyg',
                     vars = c('X', 'Y'))
## zyg: MZ
##           X1        Y1        X2        Y2
## X1 1.0000000 0.7121176 0.5930486 0.4323049
## Y1 0.7121176 1.0000000 0.5325547 0.7579110
## X2 0.5930486 0.5325547 1.0000000 0.7156957
## Y2 0.4323049 0.7579110 0.7156957 1.0000000
## -------------------------------------------- 
## zyg: DZ
##           X1        Y1        X2        Y2
## X1 1.0000000 0.7162098 0.5075549 0.3259307
## Y1 0.7162098 1.0000000 0.3067146 0.5489379
## X2 0.5075549 0.3067146 1.0000000 0.6645038
## Y2 0.3259307 0.5489379 0.6645038 1.0000000

Bivariate ACE

# Bivariate ACE model
model2 <- multivariate_ace(
  data = data2,
  zyg = 'zyg',
  ph = c('X', 'Y')
)

model2 <- mxRun(model2)

summary(model2)
## Summary of ACE 
##  
## free parameters:
##             name matrix row col      Estimate  Std.Error A
## 1  ACE.mean[1,1]   mean   1   1 -1.291566e-01 0.06079859  
## 2  ACE.mean[1,2]   mean   1   2 -1.348547e-01 0.06403822  
## 3     ACE.a[1,1]      a   X   X  1.897923e-06 0.28839284  
## 4     ACE.a[2,1]      a   Y   X  4.925751e-01 0.12721120  
## 5     ACE.a[2,2]      a   Y   Y  7.182836e-01 0.09900747  
## 6     ACE.c[1,1]      c   X   X  5.335468e-01 0.04600191  
## 7     ACE.c[2,1]      c   Y   X  2.438873e-01 0.18791282  
## 8     ACE.c[2,2]      c   Y   Y  5.127795e-01 0.13539231  
## 9     ACE.e[1,1]      e   X   X  4.191573e-01 0.02227449  
## 10    ACE.e[2,1]      e   Y   X  4.407855e-01 0.04928814  
## 11    ACE.e[2,2]      e   Y   Y  4.701544e-01 0.03325551  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:             11                    789
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:                1752.1
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/800
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:       174.0997               1774.100
## BIC:     -2428.2727               1810.381
##       |  Sample-Size Adjusted
## AIC:                 1775.504
## BIC:                 1775.532
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 04:22:46 
## Wall clock time: 0.08861446 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Output - parameters

# Output tables from the bivariate ACE
get_output_tables(model2)
## $Variance_components
##        A(%)      C(%)      E(%)
## X 0.2535936 0.3597038 0.3867026
## Y 0.5159729 0.2629640 0.2210630
## 
## $Covariation_components
##             Total         A         C         E      A(%)
## X <-> X 0.9567681 0.2426303 0.3441531 0.3699847 0.2535936
## X <-> Y 0.6861063 0.3538086 0.1250604 0.2072373 0.5156761
## Y <-> Y 0.9999194 0.5159313 0.2629428 0.2210452 0.5159729
##              C(%)      E(%)
## X <-> X 0.3597038 0.3867026
## X <-> Y 0.1822756 0.3020483
## Y <-> Y 0.2629640 0.2210630
## 
## $All_correlations
##         r_A       r_C       r_E     Total
## X <-> X   1 1.0000000 1.0000000 1.0000000
## X <-> Y   1 0.4157318 0.7246619 0.7014643
## Y <-> Y   1 1.0000000 1.0000000 1.0000000

Output - model fit

# Model fit
model2 <- ref_models(model2)

fit_stats(model2)
##   model      base ep minus2LL  df      AIC       BIC
## 1   ACE Saturated 11   1752.1 789 174.0997 -2428.273
##        CFI      TLI RMSEA  diffLL diffdf         p
## 1 1.001519 1.001072     0 16.2523     17 0.5060247

Vignette

OpenMx

A model (mxModel):

  • matrices (mxMatrix)
  • matrix algebra (mxAlgebra)
  • data (mxData)
  • fit function (mxFitFunction…)
  • expectation (mxExpectation…)
  • sub-models (mxModel)
  • CIs (mxCI), constraints (mxConstraint)

Univatiate ACE

# Twin models in OpenMx
library(OpenMx)

# Univariate ACE model
model <- mxModel(
  name = 'ACE',
  
  # Means
  mxMatrix(type = 'Full',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'mean'),
  mxAlgebra(cbind(mean, mean),
            name = 'expMeans'),
  
  # Variance components
  mxMatrix(type = 'Lower',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'a'),
  mxMatrix(type = 'Lower',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'c'),
  mxMatrix(type = 'Lower',
           nrow = 1, ncol = 1,
           free = TRUE,
           name = 'e'),
  mxAlgebra(t(a) %*% a, name = 'A'),
  mxAlgebra(t(c) %*% c, name = 'C'),
  mxAlgebra(t(e) %*% e, name = 'E'),
  
  # Expected covariance
  mxAlgebra(rbind(cbind(A + C + E, A + C    ),
                  cbind(A + C,     A + C + E)),
            name = 'expCovMZ'),
  mxAlgebra(rbind(cbind(A + C + E, .5%x%A + C),
                  cbind(.5%x%A + C, A + C + E)),
            name = 'expCovDZ'),
  
  # MZ submodel
  mxModel(name = 'MZ',
          mxData(observed = mz_data,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovMZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'X2')),
          mxFitFunctionML()),
  
  # DZ submodel
  mxModel(name = 'DZ',
          mxData(observed = dz_data,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovDZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'X2')),
          mxFitFunctionML()),
  
  mxFitFunctionMultigroup(c('MZ','DZ')),
  
  # Output tables
  mxAlgebra(A + C + E, name = 'V'),
  mxAlgebra(cbind(V, A, C, E),
            dimnames = list('X', c('Total', 'A', 'C', 'E')),
            name = 'Raw_variance'),
  mxAlgebra(cbind(A / V, C / V, E / V),
            dimnames = list('X', c('A(%)', 'C(%)', 'E(%)')),
            name = 'Variance_components')
)

model <- mxRun(model)

summary(model)
## Summary of ACE 
##  
## free parameters:
##            name matrix row col    Estimate  Std.Error A
## 1 ACE.mean[1,1]   mean   1   1 -0.01625702 0.06929625  
## 2    ACE.a[1,1]      a   1   1  0.91069985 0.09748699  
## 3    ACE.c[1,1]      c   1   1  0.37057422 0.22770404  
## 4    ACE.e[1,1]      e   1   1  0.46506705 0.03335826  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:              4                    396
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:               1067.62
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/400
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:        275.620               1075.620
## BIC:      -1030.514               1088.813
##       |  Sample-Size Adjusted
## AIC:                 1075.825
## BIC:                 1076.141
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 04:22:47 
## Wall clock time: 0.04803014 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Univariate output

# Univariate output
mxEval(Raw_variance, model)
##      Total         A         C         E
## X 1.182987 0.8293742 0.1373253 0.2162874
mxEval(Variance_components, model)
##        A(%)      C(%)      E(%)
## X 0.7010849 0.1160835 0.1828316

Bivariate ACE

# Bivariate ACE model
model2 <- mxModel(
  name = 'ACE',
  
  # Means
  mxMatrix(type = 'Full',
           nrow = 1, ncol = 2,
           free = TRUE,
           name = 'mean'),
  mxAlgebra(cbind(mean, mean),
            name = 'expMeans'),
  
  # Variance components
  mxMatrix(type = 'Lower',
           nrow = 2, ncol = 2,
           free = TRUE,
           name = 'a'),
  mxMatrix(type = 'Lower',
           nrow = 2, ncol = 2,
           free = TRUE,
           name = 'c'),
  mxMatrix(type = 'Lower',
           nrow = 2, ncol = 2,
           free = TRUE,
           name = 'e'),
  mxAlgebra(t(a) %*% a, name = 'A'),
  mxAlgebra(t(c) %*% c, name = 'C'),
  mxAlgebra(t(e) %*% e, name = 'E'),
  
  # Covariance
  mxAlgebra(rbind(cbind(A + C + E, A + C    ),
                  cbind(A + C,     A + C + E)),
            name = 'expCovMZ'),
  mxAlgebra(rbind(cbind(A + C + E, .5%x%A + C),
                  cbind(.5%x%A + C, A + C + E)),
            name = 'expCovDZ'),
  
  # MZ submodel
  mxModel(name = 'MZ',
          mxData(observed = mz_data2,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovMZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'Y1', 'X2', 'Y2')),
          mxFitFunctionML()),
  
  # DZ submodel
  mxModel(name = 'DZ',
          mxData(observed = dz_data2,
                 type = 'raw'),
          mxExpectationNormal(covariance = 'ACE.expCovDZ',
                              means = 'ACE.expMeans',
                              dimnames = c('X1', 'Y1', 'X2', 'Y2')),
          mxFitFunctionML()),
  
  mxFitFunctionMultigroup(c('MZ','DZ')),
  
  # Output tables
  mxAlgebra(A + C + E, name = 'V'),
  mxAlgebra(abs(A) + abs(C) + abs(E), name = 'Vabs'),
  mxAlgebra(abs(A) / Vabs, name = 'Aperc'),
  mxAlgebra(abs(C) / Vabs, name = 'Cperc'),
  mxAlgebra(abs(E) / Vabs, name = 'Eperc'),
  mxAlgebra(cbind(diag2vec(Aperc),
                  diag2vec(Cperc),
                  diag2vec(Eperc)),
            dimnames = list(c('X', 'Y'), c('A(%)', 'C(%)', 'E(%)')),
            name = 'Variance_components'),
  
  mxAlgebra(cbind(vech(V),
                  vech(A), vech(C), vech(E),
                  vech(Aperc), vech(Cperc), vech(Eperc)),
            dimnames = list(c('X <-> X', 'X <-> Y', 'Y <-> Y'),
                            c('Total', 'A', 'C', 'E', 'A(%)', 'C(%)', 'E(%)')),
            name = 'Covariation_components'),
  
  mxMatrix(type = "Iden",
           nrow = 2, ncol = 2,
           name = "I"),
  mxAlgebra(sqrt(I * V), name = "SD_total"),
  mxAlgebra(sqrt(I * A), name = 'SD_A'),
  mxAlgebra(sqrt(I * C), name = 'SD_C'),
  mxAlgebra(sqrt(I * E), name = 'SD_E'),
  
  mxAlgebra(solve(SD_total) %&% V, name='r_total'),
  mxAlgebra(solve(SD_A) %&% A, name = 'r_A'),
  mxAlgebra(solve(SD_C) %&% C, name = 'r_C'),
  mxAlgebra(solve(SD_E) %&% E, name = 'r_E'),
  
  mxAlgebra(cbind(vech(r_A), vech(r_C), vech(r_E),
                  vech(r_total)),
            dimnames = list(c('X <-> X', 'X <-> Y', 'Y <-> Y'),
                            c('r_A', 'r_C', 'r_E', 'Total')),
            name='All_correlations')
)

model2 <- mxRun(model2)

summary(model2)
## Summary of ACE 
##  
## free parameters:
##             name matrix row col      Estimate  Std.Error A
## 1  ACE.mean[1,1]   mean   1   1 -1.291565e-01 0.06079855  
## 2  ACE.mean[1,2]   mean   1   2 -1.348544e-01 0.06403800  
## 3     ACE.a[1,1]      a   1   1  2.990996e-07 0.28843910  
## 4     ACE.a[2,1]      a   2   1  4.925745e-01 0.12723097  
## 5     ACE.a[2,2]      a   2   2  7.182822e-01 0.09902277  
## 6     ACE.c[1,1]      c   1   1  5.335470e-01 0.04600363  
## 7     ACE.c[2,1]      c   2   1  2.438879e-01 0.18794570  
## 8     ACE.c[2,2]      c   2   2  5.127802e-01 0.13541217  
## 9     ACE.e[1,1]      e   1   1  4.191574e-01 0.02227458  
## 10    ACE.e[2,1]      e   2   1  4.407853e-01 0.04928897  
## 11    ACE.e[2,2]      e   2   2  4.701545e-01 0.03325635  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom
##        Model:             11                    789
##    Saturated:             NA                     NA
## Independence:             NA                     NA
##                |  Fit (-2lnL units)
##        Model:                1752.1
##    Saturated:                    NA
## Independence:                    NA
## Number of observations/statistics: 200/800
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty
## AIC:       174.0997               1774.100
## BIC:     -2428.2727               1810.381
##       |  Sample-Size Adjusted
## AIC:                 1775.504
## BIC:                 1775.532
## CFI: NA 
## TLI: 1   (also known as NNFI) 
## RMSEA:  0  [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 04:22:47 
## Wall clock time: 0.09123135 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.18.1 
## Need help?  See help(mxSummary)

Bivariate output

# Bivariate output
mxEval(Variance_components, model2)
##        A(%)      C(%)      E(%)
## X 0.2535930 0.3597045 0.3867025
## Y 0.5159716 0.2629651 0.2210634
mxEval(All_correlations, model2)
##         r_A       r_C       r_E     Total
## X <-> X   1 1.0000000 1.0000000 1.0000000
## X <-> Y   1 0.4157326 0.7246617 0.7014641
## Y <-> Y   1 1.0000000 1.0000000 1.0000000

Bivariate output

mxEval(Covariation_components, model2)
##             Total         A         C         E      A(%)
## X <-> X 0.9567679 0.2426296 0.3441537 0.3699846 0.2535930
## X <-> Y 0.6861056 0.3538075 0.1250609 0.2072372 0.5156750
## Y <-> Y 0.9999182 0.5159293 0.2629435 0.2210453 0.5159716
##              C(%)      E(%)
## X <-> X 0.3597045 0.3867025
## X <-> Y 0.1822765 0.3020486
## Y <-> Y 0.2629651 0.2210634

mlth.data.frame

library(mlth.data.frame)

# mlth.data.frame - multiheaded data.frame
model2 <- multivariate_ace(
  data = data2,
  zyg = 'zyg',
  ph = c('X', 'Y')
)
model2 <- def_ci(model2, ci = 'Variance_components')
model2 <- mxRun(model2, intervals = TRUE)

get_output_tables(model2)
## $Variance_components
##   A(%)-------------------------- C(%)-------------------------- E(%)-------------------------
##   est       lCI        uCI       est       lCI         uCI      est       lCI       uCI      
## X 0.2535936 0.05284689 0.5694850 0.3597038 0.070988727 0.552154 0.3867026 0.2972783 0.4942010
## Y 0.5159729 0.26094191 0.7996993 0.2629640 0.001319719 0.492202 0.2210630 0.1623230 0.3034753
## 
## $Covariation_components
##             Total         A         C         E      A(%)
## X <-> X 0.9567681 0.2426303 0.3441531 0.3699847 0.2535936
## X <-> Y 0.6861063 0.3538086 0.1250604 0.2072373 0.5156761
## Y <-> Y 0.9999194 0.5159313 0.2629428 0.2210452 0.5159729
##              C(%)      E(%)
## X <-> X 0.3597038 0.3867026
## X <-> Y 0.1822756 0.3020483
## Y <-> Y 0.2629640 0.2210630
## 
## $All_correlations
##         r_A       r_C       r_E     Total
## X <-> X   1 1.0000000 1.0000000 1.0000000
## X <-> Y   1 0.4157318 0.7246619 0.7014643
## Y <-> Y   1 1.0000000 1.0000000 1.0000000

M <- get_output_tables(model2)$Variance_components

M[['A(%)']]
##   est       lCI        uCI      
## X 0.2535936 0.05284689 0.5694850
## Y 0.5159729 0.26094191 0.7996993

# Formatting the table for html or latex report
# Powered by knitr and kableExtra
library(kableExtra)

M %>% 
  behead() %>%
  kable(digits = 3,
        align = 'c') %>%
  add_complex_header_above(M)
A(%)
C(%)
E(%)
est lCI uCI est lCI uCI est lCI uCI
X 0.254 0.053 0.569 0.360 0.071 0.552 0.387 0.297 0.494
Y 0.516 0.261 0.800 0.263 0.001 0.492 0.221 0.162 0.303 

# Saving the table for the output
M %>%
  register_output(
    name = 'Table 1',
    caption = 'Table 1: Variance components')
##   A(%)-------------------------- C(%)-------------------------- E(%)-------------------------
##   est       lCI        uCI       est       lCI         uCI      est       lCI       uCI      
## X 0.2535936 0.05284689 0.5694850 0.3597038 0.070988727 0.552154 0.3867026 0.2972783 0.4942010
## Y 0.5159729 0.26094191 0.7996993 0.2629640 0.001319719 0.492202 0.2210630 0.1623230 0.3034753

# Run at the end of the script
# to write all output tables into a single file
# openxlsx package required
write.xlsx.output(file = 'output.xlsx')

Vignettes